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A BCS (Bardeen-Cooper-Schrieffer) superconductor, which is placed out of equilibrium, can de- 
velop quantum instabilities, which manifest themselves in oscillations of the superconductor's order 
parameter (pairing amplitude A). These instabilities are manifestations of the Cooper instability. 
Inelastic collisions are essential in resolving those instabilities. Incorporating the quantum insta- 
bilities and collisions in a unified approach based on Richardson's exact solution of the pairing 
Hamiltonian, we find that a BCS superconductor may end up in a state in which the spectrum has 
more than one gap. 



INTRODUCTION 

An equilibrium BCS superconductor has the property 
that the spectral gap is located around the Fermi energy. 
When the superconductor is connected to leads thorough 
a tunnel junction and current is allowed to pass through 
the system, the super-conducting gap may be found be- 
low or above the Fermi energy. Since excitations above 
the gap are electron-like and those below the gap are 
hole-like, the shift of the gap away from the Fermi en- 
ergy is accompanied by a net difference, A A/", between 
the number of hole-like and electron-like excitations. A 
quantity 'branch imbalance' with the dimension of en- 
ergy is denoted by $ and defined by $ = where po 
is the density of states. 

Experiments [TJ [2] have directly demonstrated the dif- 
ference between the Fermi energy and the energy of the 
gap (or the condensate energy) in the same piece of super- 
conducting material. In order to understand how the in- 
jection of excitations through a tunnel junction into a 
superconductor can lead to imbalance, consider first that 
the current is injected into the superconductor as normal 
current while in the superconductor it flows as super- 
current. This means that a process by which normal 
excitations are converted into condensed pairs must take 
place in the injection region. The condensed pairs then 
flow away from the injection region as super-current. 

The processes responsible for converting normal ex- 
citations into condensed pairs are collisions, in many 
case predominately electron-phonon collisions. These 
collisions constantly convert electron-like excitations into 
condensed pairs or destroy condensed pairs by converting 
hole-like excitations. The energy of the condensate, and 
consequently of the condensed pairs, can be assigned as 
the energy of the gap (the gap midpoint). If the gap re- 
sides at the Fermi energy, particle- hole symmetry leads to 
equal rates for condensed-pair creation and annihilation, 
i.e. to no net conversion of quasi-particles into condensed 
pairs. It is the shift of the gap relative to the Fermi sur- 
face that allows for a net conversionof quasi-particles into 
condensed-pairs. 

Close to the critical temperature and under certain as- 
sumptions, to be detailed below, one can write down an 



approximate distribution function for the excitations de- 
scribing a situation where the quasi-particles are at equi- 
librium at chemical potential £i — while the condensate 
is at energy £i [3] [4] [5]: 

n (0 = : e (o + »,(o > (!) 

1 + e kT 

where e(£) and q(£) are the energy and charge of an ex- 
citation, respectively: 

^) = V / (C-6) 2 +A? ) q(0 = ^. (2) 

Ai denotes the size of the gap. Quasi- neutrality demands 

Eq. is only an approximation and has corrections, 
but for a qualitative understanding of the effect we want 
to describe, it shall be sufficient to ignore these for the 
time being. We follow here closely the treatment of [4]. 
If one plugs into the self-consistency equation, one 
obtains a relation between Ai and <£, as follows: 

A 2 = A!(0) 2 -2$ 2 , (3) 

where Ai(0) is the gap at $ = 0. It is seen that branch 
imbalance, as measured by <£, suppresses the gap (see 
Fig. (!])). In fact at $ = A(0)/V% the gap is completely 
suppressed, and all pairing correlations vanish. The sys- 
tem returns back to the normal metal state. It is known, 
however, that this state is unstable. Any normal metal 
placed below T c is unstable to Cooper pairing, i.e. to for- 
mation of pairing correlation at the Fermi energy. There 
must be a critical value of <£, designated as above 
which, the system becomes unstable to such pairing, as 
pointed out in [4]. 

The scenario associated with the instability can be de- 
scribed roughly as follows: pairing correlations are re- 
stricted to a region of size comparable to the gap, Ai 
at the gap energy, £i. As £i = becomes larger, this 
region becomes smaller because of suppression. As a 
consequence the Fermi energy remains devoid of pairing 
correlations. The cooper instability - an instability to 
the formation of correlations at the Fermi energy - then 
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conducting gap, Ai. This solution is provided by taking 
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for Ai satisfying the self consistency condition: 
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FIG. 1: Dependence of gap, Ai, on imbalance, <E>, in units 
of the gap in the absence of imbalance, Ai(0). Beyond the 
critical value of the imbalance, the systems is not stable 
for Cooper-pairing at the Fermi energy. 



takes place. 

The purpose of this letter is to analyze the state which 
is formed beyond the instability point. We shall show, 
that a novel state appears where another gap, A2, is 
formed at an energy £2 closer to the Fermi level. We shall 
make use of Richardson's exact solution of the pairing 
Hamiltonian, to describe such a state. Other approaches 
are known to be inadequate to describe the state beyond 
the instability point [6][Z]- 



SEMICLASSICAL APPROACH FOR 
INSTABILITY 



To understand better the instability that takes place, 
it is useful to consider the behavior of Anderson's pseudo- 
spins [8]: 



1 



2 I (z^°7>^/ 



h 3 



is) 



which afford a semiclassical description. The order pa- 
rameter A is given by A = A^ — iA y = | J2j s ~^i ano ^ i s a 
measure of the overall pairing correlations. The semiclas- 
sical limit of the dynamics may be obtained by recasting 
the pairing Hamiltonian as follows: 

H = J2 c L c ^> c )Ai c 3'^'^ ( 4 ) 
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through the pseudo-spins, H = X^-=i ^j s i ~ § I ^| 2 . This 
leads to the following equations of motion for the pseudo- 
spins: 



X Sn 



(5) 



A solution of d5j, where the order parameter takes the 
form A = Aier*^. Namely, the order parameter has 
a time-independent modulus, which equals the super- 



Here e(£) is given by and may be interpreted as 
the distribution function for the excitations. 

The distribution function given by ([I]) represents just 
such a solution (with time-independent modulus of the 
order parameter). The solution (|6| becomes unstable 
for <3> > The instability manifests itself in a lin- 

ear stability analysis around the solution The un- 
stable modes can be seen to represent pairing correla- 
tions forming at the Fermi energy. Namely such a mode 
contains non- vanishing components of s + (£) at energies 
around the Fermi energy. A perturbation around the 
time-independent solution will excite these modes. The 
pairing correlations around the Fermi surface then man- 
ifest themselves in a change of the order parameter A. 
Moreover the modulus of the order parameter will cease 
to be time-independent. 

If instead of ([I]) one takes the distribution function of 
the normal metal (obtained by taking Ai=0, £1 = $ = 
0) one recovers the usual Cooper instability. Such an 
instability arises if one takes a normal metal and sud- 
denly switches on the pairing interaction. This can be 
potentially realized in the lab by employing the Fesh- 
bach resonance [9[ [lOj . Such a situation was considered 
in [11]. It was shown there that the time dependence 
of the order parameter exhibits an oscillatory behavior 
which can be described as soliton trains (See figure Q). 
It should be noted that this behavior is only exhibited 
at times before collisions take place. Collisions will tend 
to equilibrate the system bringing it to the BCS state, 
with its time-independent order parameter. To treat col- 
lisions one must go beyond ([5]), however in the simple 
Cooper instability problem one already has a good quali- 
tative understanding of the system's behavior: First soli- 
ton trains will ensue, the oscillations associated with the 
soliton trains will slowly die out and the thermal equi- 
librium state, where the order parameter is constant in 
time, will be reached. 

An important development in the study of the semi- 
classical equations is the discovery of their integrability. 
This was achieved in [12] and then utilized in [13] . An in- 
teresting Wigner function formalism for the collisionless 
dynamics has appeared in [T4] . 

In the case at hand, where the Cooper instability takes 
place at the same time that a gap already exists away 
from the Fermi surface, a qualitative understanding of 
the steady state is more complicated. The semiclassical 
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FIG. 2: A sketch of the time dependence of the modulus of the 
order parameter as a function of time, the typical Richardson 
spectrum is shown for each period. After an initial quiet pe- 
riod (left Richardson spectrum), a perturbation causes a train 
of solitons to appear (middle Richardson spectrum), these os- 
cillations slowly die out and the global steady state ensues 
(right Richardson spectrum). In the Richardson spectrum 
plots, X denote instability points, while arcs are shown in 
contoured heavy lines. 



analysis just after the instability takes place was studied 
in [15]. However, without employing some new theory, no 
obvious candidate for the long time steady state beyond 
the instability point exists. Indeed does not provide 
such a solution for any value of $ and Ai. The purpose 
of the letter is to characterize this state. In particular we 
will find that two gaps coexist in the spectrum. 



QUANTITATIVE ANALYSIS OF IMBALANCE 

We shall treat the problem using a kinetic approach. 
In order for a kinetic approach to be valid, the time, 
ta = it takes for the spectrum to adjust to changes 
of the order parameter must be much smaller than the 
collision time r e . We assume a situation such as the one 
displayed in Fig. (J3|. We treat the superconductor S 
in the tunnel region. We shall later see that beyond the 
instability point two gaps appear in S as depicted in the 
figure. For simplicity Sa m &y be regarded as being in 
equilibrium. 

In addition, we shall need the assumption of local equi- 
librium, namely that the system thermalizes on a much 
faster time scale than it is forced out of equilibrium. 
The non-equilibrium forcing is associated with the injec- 
tion of imbalance into the sample. The injection process 
is balanced out in steady state by processes that relax 
imbalance, which involve excitations scattering to a re- 
gion close to the gap. We shall assume throughout that 
\T - T c \ <C T c . In this regime A«T C and only a frac- 
tion 4 of the processes involve excitations near the gap. 



This results in a typical branch imbalance relaxation rate 
which is given by Tq 1 = t" 1 ^, satisfying tq ^> r e . Lo- 
cal equilibrium establishes on a time scale r e if the non- 
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FIG. 3: Two superconductors S and Sa with Fermi energies 
Ep and Ep,A, are coupled through a tunnel junction. S devel- 
ops two gaps Ai and A 2 , at energies £1 and £2, respectively. 
Sa has gap A a- 



equilibrium forcing is small enough such that it operates 
on a time scale tq along with imbalance relaxation pro- 
cesses. 

When these conditions are satisfied namely when ta <C 
r e <C tq, one may use the following distribution function, 
n, for the occupation number of excitations, written here 
through the linear combination p = 2n — 1 for conve- 
nience: 



p(0 



+ Sp, (8) 



where e(f) = ^ - £i) 2 + A 2 , is the BCS spectrum of 
excitations and q(£) = is the charge of an excita- 

tion, and dp is the next to leading order correction in 
^. The distribution function ^ describes a situation 
where the condensate is at chemical potential £1, while 
the quasi-particles have a chemical potential — above 
it. The possibility to sustain the two entities at different 
potentials exists due to the fact that on the thermaliza- 
tion time scale r e , transitions from the condensate to the 
quasi-particle population are suppressed, since the latter 
only happen on a time scale tq. 

Quasi- neutrality dictates £1 = To find the relation 
between £1 and A one may use the self-consistency con- 
dition, J £<i£ = |, for g the interaction coupling constant 
(see 0. This gives [3] g|: 
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cT r 5T, 
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where c — 
version of 
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and 5T = T c - T. This is a corrected 
The charge imbalance may be given an 
order of magnitude estimate by balancing the injection 
rate to the relaxation rate of imbalance, Tq 1 . This is 
given by $ = ^^/r e T c 2 A _1 , for some constant b of or- 
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FIG. 4: The two branched dependence of A on the injection 
rate I, both in units of T c , for cST = -0.002, a = -1.23 and 
b = 0.276 . The state is not stable under the dashed line. 



RICHARDSON'S SOLUTION 

Richardson [16] had solved exactly the pairing Hamil- 
tonian Q The first stage for solving the Hamiltonian 
is to note that singly occupied levels do not participate 
in the dynamics and are blocked to the pairs. One may 
simply drop the singly occupied levels from the spectrum 
and solve Q for the case where levels are either occupied 
or unoccupied with pairs. 
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FIG. 5: Two dimensional plane occupied with charges. X and 
O denote unblocked an blocked levels respectively, /'s denote 
the E's. 



der unity. Here / is the injection intensity, namely, the 
inverse time it takes for an excitation to enter the su- 
perconductor. The contribution of the integral in ([9| 
can also be estimated based on a expansion of the ki- 
netic equation by order of ^ The result is aT c 2 /r e , 
for some constant a of order unity. This constant may 
be either positive or negative depending on whether the 
injection process tends to suppress or enhance [2 super- 
conductivity, respectively. We shall assume a situation 
where a is negative [4 j, namely enhancement [20] . 

Putting the above estimates into (|9| one obtains a re- 
lation between the gap A and the injection rate: 



2A 2 
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{bl)\ A 2 =cT c 5T-aI, (10) 



where / = IT^r e . For a given / there exist two different 
A corresponding to the + and — signs in ( [To] ) . Using a 
stability analysis of the semiclassical equations ([5| it was 
shown in [4 that the superconductor becomes unstable to 

A 2 

cooper pairing for A < -^ E , this occurs when the injection 
rate is above some critical value, and only for the solu- 



tions in (10) corresponding to the minus sign. It was also 



shown in [4] that beyond this point the order parameter 
starts oscillating on the time scale ta- This invalidates 
the distribution function ([sj) along with the kinetic ap- 
proach based on the BCS expressions for the spectrum 
e(£) and excitation charge <?(£)• The BCS approach is 
no longer valid as it is based on the assumption that the 
order parameter is time independent on a time scale ta- 
We shall argue in the sequel that ([8| is in fact valid, al- 
beit with a new spectrum of excitations, consisting of two 
gaps, and a suitable both computed making use of 
Richardson's exact solution of the problem. 



We take the Hamiltonian Q, where j runs over 2M 
levels around the Fermi energy. Suppose that there are 
N + 2M particles in the system and L levels are blocked. 
Richardson's solution for the eigenstates for the reduced 
spectrum is the following: 

N + 2M-L 

fl 62.10), bi = j2^iAA,i (ii) 



Where the E"s are comp lex parameters which must satisfy 
the following non-linear equations: 



E v — En' 



(12) 



which can be viewed as the condition of vanishing elec- 
tric field (equilibrium) for each of the charges E u , in two- 
dimensions. Each E carries a charge +1, while at each 
location of & a charge — ^ is present. In addition, there 
is a background field, — ^. An equilibrium point where 
E can be placed can always be found between any two 
unblocked £'s. In the continuum limit, where the level 
spacing is the smallest energy scale, this defines a charge 
density, A(£), on the real axis, which is composed of the 
contribution of both the i?'s and the unblocked £'s. As- 
suming a constant density of state po for the <^'s, we may 
also define = 2p^ 1 A(^). In addition to the E's on 
the real axis some of the £?'s arrange themselves in the 
complex plane. In the continuum limit the distribution 
of the complex E's is given by line densities on arcs in 
the complex plane [17 [18 . The distribution of the E's 
may be termed as the 'Richardson spectrum' (See Fig. 
(|5|). The position and shape of the arcs is determined by 
electrostatic equilibrium. may be understood as the 
generalization of the occupation number of excitations in 
the BCS approach. Indeed, To effect an excitation one 
must change the density 
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For a given and assuming m arcs in the complex 
plane harboring the complex £"s, a continuum limit so- 
lution for the charge distribution is given in the following 
form H3 HH]: 



HO 



\ 



The jump discontinuity of the electric field are caused by 
the line density of the charges. The jump discontinuities 
of h consists of a jump discontinuity on the real axis, and 
jump discontinuities on the branch cuts of R(£) which 
are drawn as curved lines stretching from & + ^A^ to 
£i — iAi. The values of ^ and A^ are constrained by 
self-consistency conditions: 



/ 



gg = 2, 

R(0 9 C 



Km -I. 



(13) 



Note that for m = 1 we get the usual self-consistency 
condition, where p(£) is to be identified with 2n(£) — 1. 

The energy of the Richardson state is given by ^2 U E v + 
Y^j iij, where £ tij are the blocked states. As h(£) encodes 
the density of £"s, it is possible to compute the energy 
and the number of particles (here = R(C)) : 



E = 



I 



R(0 



I 



[R(Q] + p(0 
R(0 



di (14) 



where [...]+ denotes taking the non- negative (polyno- 
mial in £) part of the Laurent expansion around infinity. 



QUASI-PARTICLES AROUND INSTABILITY 
POINT 

To describe the Cooper instability in the framework of 
Richardson's approach. Consider a normal metal state 
placed below T c . Since A = no arcs exist and all the E"s 
are on the real axis, distributed such that the total charge 
density corresponds to Fermi distribution. The instabil- 
ity manifests itself as a vanishing of h(£) at two com- 
plex conjugated points. Since the field vanishes, these 
points may be occupied by £"s. Indeed, collisions will 
cause these to start populating until a full arc is formed, 
describing the equilibrium superconducting state, as de- 
scribed in Fig. (2). 

By inserting (8) into the the expression for /&(£), one 

A 2 

sees that points of vanishing field appear for A < -rf- 
(this conclusion was reached by semiclassical means in 
[4]). The resolution of the instability is achieved when 
an arc appears near the Fermi energy, in addition, that 
is, to the arc at £i. 

We now study the state with two arcs, m = 2. By 
Hellmann-Feynman, q(£i) = f|s assuming that the dis- 
tance between the two arcs is much larger than their 



size, one gets: q(£) = — ^ . Note that this form 

of q(£) dictates the following excitation representation: 
Both well below £i and well above £2, = 1 an d ex- 
citations are electron-like. In the region between £1 and 
£2 but well away from both gaps, q(£) = — 1, and excita- 
tions are hole-like. This choice does not agree with the 
standard choice of the excitation representation. How- 
ever, the choice of excitation charge far away from the 
gaps is a matter of representation. 

Consider now transitions due to collisions between 
states characterized by two arcs. To compute the tran- 
sition rates, one must compute the matrix element of 
the phonon interaction between two different Richardson 
states. We are interested in solving the kinetic equation 
only to first order in ^, so one may assume that one of 
the excitations involved in the transition, has energy of 
order T, and thus has the same character as in the nor- 
mal metal. Under such an assumption, the transition of 
this excitation to a level near the gap is dictated by the 
one particle density matrix of the target level. The one 
particle density matrix consistent with computed 
above, can be easily found. From its form, it is easy 
to see that, the transition rates may be defined by coher- 
ence factors just in the BCS case, and these are given by 

The other ingredient in a kinetic approach is the energy 
conservation delta function which features in the Fermi 
golden rule. To study energy conservation we make use of 
the following picture valid close to T c due to Pethick and 
Smith [7J, which states that a change in the number of 
particles may be divided into a contribution of the change 
of the normal charge, Sq n = f qSpd^ and superconduct- 
ing charge, Sq s = j p5qd£, as follows SN = Sq s + 5q n . 
Starting from the expression for the energy of the level, 
we now write down the variation with respect to a change 
in the density p of the combination W = E — vN. Con- 
sider the energy of an excitation next to ^: 



SW 



yU-£) 2 + A^ 



the (— Y arises from the unusual choice of the excitation 
representation discussed above (we assume £1 < £2 < 0). 
The energy of an excitation, neglecting the charge trans- 
fer to the condensate, is written as e u (^) = 

P 5q s =0 

The kinetic equation may now be written and solved 
taking into account the energy change due to charge 
transferred to the condensate (Sq s ). It turns out how- 
ever, that these contribute to a lower order. As a re- 
sult the distribution function turns out to be = 

tanh ^ e ^2T^ ) + 3p ^ s f° un d by balancing the injec- 
tion rate with the relaxation provided by the linearized 
collision operator. The linearization is performed around 
the normal-metal state, so the relaxation times are inde- 
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FIG. 6: The dependence of the two gaps on the injection I. 
The dashed lines are £i and £2 (£2 lies above £1). The solid 
lines are (from bottom to top) , £1 — Ai , £1 + Ai , £2 — A2 , £2 + 
A2, and delineate the gaps. The approximations leading to 
this solution become poorer as the gaps approach one another. 



pendent of A, which yields that dp in (JSJ) depends only 
on / and temperature. 

To obtain £1 and £2 we separate the contribution of the 
super-current from the region close to £1 in energy and 
the region close to £2 . We demand that the charge carried 
away as supercurrent associated with a particular region 
is equal to the amount of charge transferred to the con- 
densate by collision into that region. The superconductor 
velocity v s is the same regardless of the region considered, 
while the contribution to the density of superconducting 
pairs from region i is proportional to Af. Collisions to 
region i transfer charge to the condensate proportional 
to £^A^ . As a result, one obtains £^ 



2^2 A^+A| 



Plugging these estimates into the self consistency equa- 



predicted. The multi-gapped state also manifests itself 
in oscillations of the order parameters. Experiments in- 
volving the Josephson junction coupled to the injection 
region will be sensitive to these oscillations. The form 
of oscillations can be found using the quantum to classi- 
cal correspondence between Richardson's states and the 
semiclassical solutions found in [T2j . Applying this rela- 
tion to an exact description of the oscillatory behavior 
will be the subject of future work, however without fur- 
ther analysis one can conclude that, due to the fact that 
the system describes two condensates at energies sepa- 
rated by a distance ££, the oscillations will predominately 
be at a frequency of We have developed the theory 
in the simplest case where the two condensates are well 
separated at a temperature near T c . It is interesting to 
extend the theory beyond this simple application. This 
may be done as a kinetic theory is possible due to the fact 
that the transition rates between Richardson states and 
the energy of any such state is known or feasibly com- 
putable [19]. The complications of solving the kinetic 
theory may be overcome by using a numerical approach. 

In the specific situation discussed in this paper two 
spectral gaps appeared. The theory, however, extends to 
cases where more than two gaps exist. Roughly speaking, 
additional gaps appear at points, £0, where n(£o) = 1/2 
and pairing correlations are small enough as not to sup- 
press the appearance of a new gap. Experimental setups 
can potentially be created that would deform the distri- 
bution function so strongly as to create such points (the 
distribution function must be made to pass through 1/2 
and the other gaps must be suppressed). Given the dis- 
tribution function, a quantitative treatment can be given 



using the generalized self-consistency conditions, (13). 



tions (13) one obtains equations relating the gaps to the 



current intensity: 
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(15) 



where 1 = 2 and 2 = 1. A sketch of the results for the 
same parameters for which Fig. Q is drawn is given in 

Fig. del). 
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CONCLUSION 

We have seen that non-equilibrium effects may excite 
a mode by which a superconductor develops another gap. 
The possibility of having two gaps is most easily revealed 
by considering the multi-arc Richardson states. These 
gaps will manifest themselves in the same type of experi- 
ments that reveal the BCS gap. A probe for the density of 
states (such as another tunnel junction) must be coupled 
to the injection region. In situations where the current 
exceeds the instability point, a multi-gapped structure is 
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